#include "base.h"


int main(int argc,char **argv)
{


          //////////////////

          char kappa_sea[1024];
          int Nmass;
          char input[1024];
          char output[1024];
          char P_to_use[1024];
          int  set_p[4];
	  int iprop;
          int twisted_basis;

          int beta100;

          char momentum_order[1024];
          char gamma[1024];
	  int r;

          ////// Parsing command line

          for(int i=1; i<argc; i++) {
            if(strncmp(argv[i],"-kappa_sea",10)==0)  sscanf(argv[i],"-kappa_sea=%s",kappa_sea);

            if(strncmp(argv[i],"-Nmass",6)==0)  sscanf(argv[i],"-Nmass=%d",&Nmass);
            if(strncmp(argv[i],"-input",6)==0)  sscanf(argv[i],"-input=%s",input);
            if(strncmp(argv[i],"-output",7)==0)  sscanf(argv[i],"-output=%s",output);
            if(strncmp(argv[i],"-T",2)==0)  sscanf(argv[i],"-T=%d",&T);
            if(strncmp(argv[i],"-L",2)==0)  sscanf(argv[i],"-L=%d",&L);
	    if(strncmp(argv[i],"-r",2)==0)  sscanf(argv[i],"-r=%d",&r);
            if(strncmp(argv[i],"-beta100",8)==0)  sscanf(argv[i],"-beta100=%d",&beta100);
            if(strncmp(argv[i],"-P_to_use",9)==0)  sscanf(argv[i],"-P_to_use=%s",P_to_use);
            if(strncmp(argv[i],"-iprop",6)==0)  sscanf(argv[i],"-iprop=%d",&iprop);
            if(strncmp(argv[i],"-set_p",6)==0)  sscanf(argv[i],"-set_p=%d,%d,%d,%d",&set_p[3], &set_p[2], &set_p[1], &set_p[0]);


             if(strncmp(argv[i],"-momentum_order",15)==0)  sscanf(argv[i],"-momentum_order=%s",momentum_order);
             if(strncmp(argv[i],"-gamma",6)==0)  sscanf(argv[i],"-gamma=%s",gamma);

             if(strncmp(argv[i],"-twisted_basis",13)==0)  sscanf(argv[i],"-twisted_basis=%d",&twisted_basis);

          }
          int Nmass_comb=Nmass;
          for (int i=1; i<Nmass; i++)  Nmass_comb=Nmass_comb+(Nmass-i);
          init_base_gamma(gamma);

          printf("==============================\n");
          printf("Nmass=%d\n",Nmass);
          printf("Nmass_comb=%d\n",Nmass_comb);
          printf("L=%d T=%d\n",L,T);
          printf("set=(%d,%d,%d,%d)\n",set_p[0], set_p[1], set_p[2], set_p[3]);
          if (twisted_basis) printf(" Propagators in base twisted\n");
          else  printf(" Propagators in physical base \n");
          printf("input=%s \n",input);
          printf("output=%s \n",output);

          printf("==============================\n");
  

	  char momentafile[1024]; sprintf(momentafile,"mom_%d.dat",beta100);
	  read_momenta_data(momentafile);
	 double *nn=(double*)malloc(sizeof(double)*(Nmom));
	  compute_nn(momentum_order,nn);
	  Momentum_Magro(momentum_order,set_p);


	  double *P2tilde=(double*)malloc(sizeof(double)*(Nmom));
	  double **Ptilde=(double**)malloc(sizeof(double*)*(Nmom));
	  for (int imom=0 ; imom<Nmom; imom++) Ptilde[imom]=(double*)malloc(sizeof(double)*(4));
//	  for (int imom=0 ; imom<Nmom; imom++) printf("%8.5f\n",SinP2[imom]);
	  for (int imom=0 ; imom<Nmom; imom++) P2tilde[imom]=SinP2[imom];
	  for (int imom=0 ; imom<Nmom; imom++)  for (int idir=0; idir<4; idir++) Ptilde[imom][idir]=SinP[imom][idir];
	  if(strcmp(P_to_use,"P")==0)
  	{
                 printf("Using P insteadd SinP\n");
                  for (int imom=0 ; imom<Nmom; imom++) P2tilde[imom]=P2[imom];
                  for (int imom=0 ; imom<Nmom; imom++)  for (int idir=0; idir<4; idir++) Ptilde[imom][idir]=P[imom][idir];

  	}


         
	  colorcolorspinspin **s0FT_r=malloc(sizeof(colorcolorspinspin*)*(Nmass));
          for(int imass=0; imass<Nmass; imass++) s0FT_r[imass]=malloc(sizeof(colorcolorspinspin)*(Nmom));

         char s0ft_input[1024];
	 sprintf(s0ft_input,"%s",input);
	 printf("%s\n",input);
         FILE *fin=fopen(s0ft_input,"rb+");
         if(fin == (FILE *)NULL) {
                        fprintf(stderr,"Unable to open file %s for reading\n",s0ft_input);
                        return EXIT_FAILURE;
            }
         for(int imass=0; imass<Nmass; imass++)
                   {
                  // printf("im=%d\n",imass);
                   for (int imom=0; imom<Nmom; imom++)
                        for (int ic_ext=0; ic_ext<3; ic_ext++) for (int is_ext=0; is_ext<4; is_ext++) for (int ic_int=0; ic_int<3; ic_int++) for (int is_int=0; is_int<4; is_int++){
                        fread(&s0FT_r[imass][imom][ic_int][ic_ext][is_int][is_ext],sizeof(complex),1,fin);
                        //printf("%8.10f +i %8.10f\n", creal(s0FT_r[imass][imom][r][ic_int][ic_ext][is_int][is_ext]),cimag(s0FT_r[imass][imom][r][ic_int][ic_ext][is_int][is_ext]));
                        }
                   }
           fclose(fin);


	printf("Read\n");


	  //form factors
 	 complex ***s0_ff=malloc(sizeof(complex**)*(Nmass));
 	 for (int imass=0; imass<Nmass; imass++) s0_ff[imass]=malloc(sizeof(complex*)*(Nmom));
	 for (int imass=0; imass<Nmass; imass++) for (int imom=0; imom<Nmom;imom++) s0_ff[imass][imom]=malloc(sizeof(complex)*(3));

	for (int imass=0; imass<Nmass;imass++)
	{
	  for (int imom=0; imom<Nmom; imom++) 
          compute_s0tw_ff(s0FT_r[imass][imom],s0_ff[imass][imom],Ptilde[imom],P2tilde[imom],nn[imom],r);


	  char s0ff_output[1024];
	  FILE *fout=NULL;
	  for (int sigma=0; sigma<3; sigma++)
 	 {

        	// ff S(p)

	         sprintf(s0ff_output,"%s_im%d_sigma%d.out",output,imass,sigma+1);
        	 fout=fopen(s0ff_output,"w");
	          if(fout == (FILE *)NULL) {
        	  fprintf(stderr,"Unable to open file %s for writing\n",s0ff_output );
	          return EXIT_FAILURE;
    		  }

	  for (int imom=0; imom<Nmom; imom++)
	     fprintf(fout,"%8.5f %8.10f + I %8.10f\n",P2tilde[imom],creal(s0_ff[imass][imom][sigma]), cimag(s0_ff[imass][imom][sigma]));
	   fclose(fout);
        }

       }

  return 0;
}
